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ABSTRACT 

Instabilities of the dust layer in a protoplanetary disk are investigated. It is known 
that the streaming instability develops and dust density concentration occurs in a situ- 
ation where the initial dust density is uniform. This work considers the effect of initial 
dust density gradient vertical to the midplane. Dust and gas are treated as different 
fluids. Pressure of dust fluid is assumed to be zero. The gas friction time is assumed to 
be constant. Axisymmetric two-dimensional numerical simulation was performed using 
the spectral method. We found that an instability develops with a growth rate on the 
order of the Keplerian angular velocity even if the gas friction time multiplied by the 
Keplerian angular velocity is as small as 0.001. 

This instability is powered by two sources: (1) the vertical shear of the azimuthal 
velocity, and (2) the relative motion of dust and gas coupled with the dust density 
fluctuation due to advection. This instability diffuses dust by turbulent advection and 
the maximum dust density decreases. This means that the dust concentration by the 
streaming instability which is seen in the case of a uniform initial dust density becomes 
ineffective as dust density gradient increases by the dust settling toward the midplane. 



Subject headings: planetary systems: protoplanetary disks- 
hydrodynamics — instabilities 



-solar system: formation- 



1. INTRODUCTION 

The first step for planetary formation requires the formation of planetesimals larger than km- 
size. However, the mechanism that planetesimals are formed from dust is scarcely understood. 
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one is the continuous growth to planetesimals due to t he sticking (jWeidenschilhng &: Cuzzil 119931 : 
Cuzzi et al.lll993l : IWurm et al-lboOll : kornet et allboOlh . 



The gravitational instabihty has advantage because the problem that meter-sized dust falls 
toward a central star is avoided. The gravitational instability is theoretically more tractable be- 
cause the formation process dose not depend on poorly-understood surface forces of dust material. 
However, if the disk is turbulent, the critical density for the gravitational instability cannot be 
reached because the dust is stirred from the midplane. The magneto rotational instability (M RI) is 
a candidate for turbulent sources in the stage of the planet formation (jBalbus &: Hawlevlll99ll ) . The 
ionization degre e is so low in the planet forming region that there is possibility that the MRI does 
not occur there (jSano et al.ll200Cll ). It is pointed that the dust layer can become t urbulent because 



of th e shear instability even if the global turbulence such as MRI do es not occur (IWeidenschilline 



1980 ) . The shear instabili ty has be en ex t ensively studied analytically (ISekiyalll998l : ISekiYa &: Ishitsu 
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Recently, a different type of instability in the dust layer draws attention. The gas supported 
by a negative radial pressure gradient revolves slower than the Kepler velocity. The dust revolves 
against the head wind. As a result, the dust falls tow ard the central star by losing the angular 



momentum (lAdachi et al 



1976 



WeidenschillingI 119771'). On the other han d, the gas moves out 



ward in the disk by gaining the angular momentum (jNakagawa et al.lll986l ). lYoudin &: Goodman 
(|2005l ) has fo und the streaming instabi lity c an occur when dust and g as have relative velocities. 
Furthermore. lYoudin &: JohansenI (j2007l ) and I Johansen &: YoudinI (j2007l ) have detailedly performed 
analysis and numerical simulations on the streaming instability in the situation that the unper- 
turbed dust density is spatially homogeneous and the gravity in the axial direction of rotation is 
ignor ed. The simulations sh owed that the streaming instability has the action which concentrates 
dust. iJohansen et al.l (j2007l ) has presented that Ceres-sized planetesimals are formed by the dust 
concentration due to the streaming instability and the self-gravity when there are meter-sized boul- 
ders in the MRI turbulence disk. However, the concentration by the streaming instability is only 
effective for dust larger than 10cm in size. It is not known whether dust can grow up to this size. 

It has been yet to be studied what instability occurs when the dust has distribution with 
a vertical gradient. It is important to examine instability in this case in order to understand 
the growth of dust. In this work, we perform two-fluid of gas and dust with a single size, two- 
dimensional simulations. We present that instability occurs and flow transits into turbulence if 
dust with cm-sized has a graded density distribution. 

In §2, the formulation is performed under assumptions of the fluid approximation of gas and 
dust. In §3, numerical results are presented. In §4, we discuss the energy sources of instability by 
deriving energy equations from linearized equations. In §5, we conclude. 
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2. BASIC EQUATIONS 



This section gives the basic equations. The gas friction force is characterized by the friction 
time Tf, which is the time during which the relative velocity of a dust aggregate and the gas becomes 
1/e. The friction time depends on the radius of dust aggregate a, the dust solid density ps, the gas 
dens ity Pg, and th e thermal velocity c^h- In the Hayashi mod el, the friction time is given Epstein's 
law (|Epsteinlll924l ) and Stokes' law (ILandau fc Lifshitglll987l ) 
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where p is the mean molecular weight, rriH is the mass of a hydrogen atom, and o"moi is t he mean 
cross section of the molecules, fg is the gas density ratio compared to the Hayashi model (jHavashi 
198lh . 



We assume that all dust aggregates have an identical friction time, and treat dust aggregate s 
as a pressure-less fluid. The latter assumption is good only if tjQx 1 (|Garaud &: LinI 12004 ). 
However we performed numerical simulations for wide rage of the value tj^k in order to understand 
basic physics of dust-gas two fluids. 

We neglect the curvature of the cylindrical coordinates (r, 0, z) and use the local Cartesian 
coordinate system which rotates with the Keplerian velocity. Our coordinates x, y, and z denote 
radial, azimuthal, and vertical directions of the disk, respectively. That is, x = r — R, y = 
R[(j) — and z, where Qk^R) is the Keplerian angular velocity at a fiducial radius r = R 

and we neglect higher order terms of x, y and z. In the following, we denote vk{R) and ^k{R) 
by vk and VLk for simplicity. The gas can be assumed to be incompressible because the dust 
layer treated here is much thinner than vertical scale height of the gas disk and, in addition, the 
flow velocity is subsonic. The vertical components of the gravity of the cent ral star and the disk 
self-gravity are neglected; this assumption is also used in t he previous works (lYoudin k. Goodman 
20051 : lYoudin &: JohansenI 120071 : iJohansen &: YoudinI 120071 ) . We consider an unperturbed state in 
which the radial pressure gradient dPo/dR is a negative constant. 

We assume axisymmetric flows, i.e. physical quantities are independent of y. Thus, we obtain 
the continuity equations, the momentum equations of gas and dust, 
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^ + V2 • {pdUd) = yoVW (6) 

-K^ + {Ud- V2)Ud = 2Ud X + 'i^lxx --{Ud- Ug), (7) 

where V2 = (^;0, fix = (0, 0,i7i<-), P is the local gas pressure perturbation, pd is the 
dust density defined by the total dust mass floating in a unit volume, and U = {U, V, W) is 
velocity of gas and dust, the subscripts g and d denote gas and dust, respectively. In order to 
solve the continuity equation of dust stably, the diffusive term is added to equation([6]) artificially. 
The diffusive parameter is chosen so that the numerical stability is maintained, and also the 
numerical diffusion is negligibly small. We confirmed that dust mass was conserved in this method. 



We perform velocity translation in the azimuthal direction following iJohansen et al.l (l2006al ). 
The system which balances between global pressure gradient, centrifugal force, and gravity at r = i? 
rotates with a sub-Kepler velocity given by 

3 

Vq = --^KX - r]VK, (8) 

where 77 is a dimensionless parameter which expresses the effect of global radial pressure gradient 
and defined by 

1 f)P 

If the reaction force of dust on gas through the friction is negligibly small, the gas revolves at tjvk = 
54 m s~^ slower than the Kepler velocity in the Hayashi model disk. Substituting U = u — Voy 
into equations dH) - ([7|) gives 

V2-Ug = 0, (10) 

^ + {Ug ■ V2)Ug = -— V2P + 2Ug xnK + ^Ug^KV " {u g " Urf), (11) 

dt Pg 1 TfPg 

^ + V2 • {pdUd) = vdVW (12) 
+ {Ud ■ V2)Ud = 2Ud xftK + ^UdO.Ky - —{ud - Ug) - 2t]Vk^kX. (13) 
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Equations (jl0p - (jl3p are solved with the Fourier spectral method. Boundary conditions are periodic. 
As for periodic boundary condition of the z direction, computational box size is large enough 
for dust not to cross the boundary and for the eigenfunction obtained by the normal mode analysis 
to decay. A second-order Adams-Bashforth scheme for the non-linear terms and a Crank-Nicolson 
schem e for the viscous te rms are employed. We use the phase shift method to eliminate aliasing 



error (jCanuto et al.lll988l ). 
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2.1. INITIAL CONDITIONS 



Youdin fc GoodmanI (|2005l ) , lYoudin JohansenI (|2007l ) and Ijohansen fc YoudinI (|2007l ) inves- 
tigated the streaming instability of the dust layer in a protoplanetary disk using an unstratified 
uniform dust density distribution as their initial conditions. However, a dust concentrated region 
in which the dust density has the same orders of magnitude with the gas density would be realized 
if dust settle toward the midplane, and the dust density should decrease with the distance from the 
midplane. As a simple model of this setting, we employ the initial density distribution which has 
constant dust density around the midplane and sinusoidal transition zones: 
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where Zd the half-thickness of the dust layer, and hd the half-thickness of the transition zones, 
where the dust density varies from PdoiO) to sinusoidally. Here the half-thickness of the dust layer 
is given by 



2pdi0) 

and the surface density of the dust is given by 
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7.1/d(r/AU)"i-^ g cm^ for r < 2.8AU, 
30/d(r/AU)"i-^ g cm^ for r > 2.8AU, 
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where fd is a parameter (fd = 1 f or the Hayashi model). We used Hayashi's solar nebula model 
(|Havashilll98ll : lHavashi et al.lll985l l at lAU as the d ust surface density S^. Initial velocities of dust 
and gas are given quasi-stationary flow obtained by iNakagawa et al.l (|l986l ) . 
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Note that the co-ordinate system in this paper mov es with sub-Kepler velo city (1 — r])vK', on the 
other hand, it moves with the Kepler velocity vk in INakagawa et al.l (|l986l ). 
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3. NUMERICAL RESULTS 

This section displays results of the numerical simulation of dust and gas, two-fluid simulations. 
Model parameters used in this work are listed in Table 1. 



3.1. Constant Dust Density Distribution 



We now present a simulation result under the constant dust distribution condition for com- 
paring our work to lYoudin JohansenI (j2007l ) and iJohansen &: YoudinI (|2007l ) . Figure [1] shows 
snapshots for the evolution of dust density, where Pd (0)/pg = l,hd/zd = 0, and tjQk = 1- The 
streaming instability occurs similar to the results of iJohansen &: YoudinI (120071 ) . The dust con- 
centration is seen after turbulence fully develops. Table 1 shows the growth rate and the 
radial wave number k of the most unstable mode. We performed the linear analysis based on 



Youdin &: JohansenI (j2007l ). and confirmed that the streaming instability is reproduced well with 
an accuracy of 1 % error in the growth rate. 



3.2. Stratified Dust Density Distribution 



Here, results for stratified dust density distributions are presented. First, we investigate the 
case of TfClx = 10~^. This friction time corresponds to approximately 1cm dust diameter at 1 AU 
from equation ([T]) . Figure [2] shows the initial profiles of the dust density and velocities of gas and 
dust for pd{0)/pg = 1 and ha/za = 0.5. It is seen that the radial velocities of dust and gas are very 
slow, \ug\/{r]Vk), \ud\/{i]Vk) ~ 10~^. Figure [3] shows snapshots of the evolution of dust density. The 
perturbation of short wave length grows by an instability, and the flow transits to turbulent state. 
And then, the dust is stirred from the midplane. Because the growth time of instability ~ lyr 
is much shorter than the settling time of dust l/^TjQ'j^) ^ 10^ yr, the result would not change 
even if there is the vertical gravity which is omitted in our simulation. Thus, this instability may 
prevent the disk from the gravitational instability. 

The dust concentration occurs as a result of "strea r ning i n stability" if the dust density dis- 
tribution is constan t as s hown by lYoudin Goodman! (|2005l ). lYoudin &: JohansenI (j2007l ). and 
Johansen YoudinI (j2007l ). However, if an initial dust density is not homogeneous, maximum dust 
density does not exceed initial one. The growth rate of instability is approximately the Keplerian 
frequency as seen in Table 1. The detail of the instability in this case is analyzed in §4. Figure H] 
shows the evolution of r.m.s. of vertical dust velocities weight-averaged by dust density {wd)rms- 
Here 
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(several m/sec). 

Next we consider the case of tjQk = 1- This friction time approximately corresponds to dust 
with a = Im at lAU from equation ([2]). In this case, we expect that the random velocities of 
dust is large, so that our approximation of pressure-less fluid for dust is not justified. However, we 
perform calculations for this case in order to understand the essence of the dust and gas two-fluid 
instability by comparing the results with that for tjQk = 10~^. Figure [5] shows the initial profiles 
of dust density and velocities of the gas and dust for P(i{0)/pg = 1 and h^/zd = 0.5. A m-sized 
body falls toward a central star with velocity r]VK- As the dust settles toward the midplane to form 
the dust layer, dust density grows compared with the gas density. Then, the infall velocity of dust 
toward the central star slows down due to the increase of dust inertia as seen from equation (jl9p . 
On the other hand, gas moves outward by obtaining the angular momentum from dust as seen from 
equation (fTTj) . Radial velocities of dust and gas vary with the distance from midplane depending 
on the dust density; that is, the vertical shear of radial velocities arises. The density pattern of the 
shear instability is actually seen at the upper right panel in Figure [5l It is shown in §4 that the 
energy source of this instability is vertical shear of the azimuthal and radial components of the gas 
velocity. 

The r.m.s. of a vertical dust velocity for tjQk = 1 is two orders magnitude smaller than that 
for Tfilx = 10~^ as se en in Figure HI The r e ason would be that the coupling between gas and dust 



is weak. As shown by iJohansen &: YoudinI (j2007l ). dust concentrates by the streaming instability 



after the flow becomes turbulent. That is also seen in our simulation at tVtK = 26.5 as shown in 
Figure [H 

In our calculation, the dust does not settle down because the gravity of central star is neglected. 
In reality, the dust would settle down before the start of the shear instability for the parameter, 
Pd{0)/pg = l,hd/zd = 0.5 and tjUk = 1 because the settling time l/(rjO|^) ~ is shorter 

than the growth time of the turbulence. As dust settling proceeds, the shear instabil ity occurs when 



the growth time of shear instabi lity becomes shorter than the dust settling time (jJohansen et al 



20n6bl : l.lohansen fc YoudinI 120071 ^ 



4. ENERGY SOURCES OF INSTABILITY 

In this section, we analyze the energy source of the instabilities presented in §3. The energy 
budget for the instability can be estimated from distributions of density and velocities obtained by 
our simulations shown in §3. 
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4.1. Linearization 

We linearize equations p^ - (fT3|) as follows: 

du' dw' , ^ 



5i 9x ^ dz Pg dx " TfPg " TfPg 

dv' _ dv' dVg 1 ^ , 1 - / / l^ 1 



+ %^ + K = -2^^^9 - T^^^^''^ - ^rf) - 7::r^d(^s - ^d), (25) 

Assuming /' tx q^^^^-^^) [k denotes the radial wave number, and uj = lor + iuji denotes the complex 
frequency), equations ([23|) -([30 |) are written as 
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- iCjdu'd + w'd^ = ^^Kv'd -^{u'd- Ug), (36) 
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where ujg = lo — kug and uid = to — ku. 
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4.2. Energy Equations 



We consider the radial energy budget for gas. Multiplying equation (j32|) by u * yields 
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where the superscript * denotes the complex conjugate. Taking the real part yields 
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We perform the similar manipulations for azimuthal and vertical directions, 
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Similarly, we consider the energy budgets for dust. 
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Here, let us see the total energy budget. Multiplying equations (|l0l) - (P2]) by pg and equations 
(P3|) - (P5]) by pdi we sum them. Note that equations (jH]), (HH) is multiplied by 4 in order to eliminate 
Epicyclic terms. 
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The equation (j46p is integrated in the z-direction in order to examine the total energy budget. We 
focus on the pressure term of the second term of the right hand in the equation (j46p . 



— — w„ dz 

dz ^ 



dz 



dz 



P'iku'gdz, 



(47) 



where the partial integration and the periodic boundary condition are used in the first equality, 
and the equation (i3T]) is used in the second equality. We have 



P'iku'* + 



dP' 



dz = 0. 



(48) 



Thus, the total energy budget of the system is determined by the shear terms, the friction terms 
due to the fluctuation of the dust density, and the frictional dissipation terms; the total energy does 
not change by the pressure terms and epicyclic terms, and they merely exchange energy among x, 
y, and z directions. Figure [7| shows the each term in energy equations (I40p - (145p and Figure [8] shows 
the each term in total energy equation ()46p for krjr = 10.8 at which the growth rate has maximum 
value for TfQx = 1 and pd{0)/pg = 1. Then the shear of radial and azimuthal velocities of gas is a 
main energy source. 



Figure [9] shows the each term in energy equations (|40 p - (j45p and Figure [10] shows the each 
term in total energy equation (j46[) for krjr = 48.6 at which the growth rate has maximum value for 
TfQ.K = 10~^ and pd{0)/pg = 1. The shears of azimuthal velocities dvg/dz and dvd/dz are main 
energy sources. In addition, the radial relative velocity between dust and gas, coupled with the 
density fluctuations, are also an important energy source. 

In energy equation ([l6|) . the last term {{^g — Ud)^ [p'd'^'g] + ^i^g ~ ^d)^ [Pd''^'g] } shows 
an instability powered by relative velocity of gas and dust, coupled with Eulerian dust density 
fluctuation. Substituting the linearized continuity equation for dust ([35]) . 
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this term is divided into parts 
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(50) 



The first term in the right hand of equation (I50p presumably corresponds to th e streaming instabil- 
ity addressed by lYoudin &: GoodmanI ([20051 ) , lYoudin &: JohansenI ([20071 ) and [Johansen &: Youdin 

dw' \ 

iku'^ + I denotes that the increase of the fluc- 
tuation of dust density is related to the streaming instability. If pg ^ pd and Tfilx ^ 1) the 
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divergence of the dust velocity approaches zero owing to the incompressibihty of gas and a small 
relative velocity between dust and gas. Thus, if dust size is small, the contribution of the first term 
of the rig ht hand side of equat i on (1501) to the instability becornes srn all. This is consistent with the 
results of lYoudin &: GoodmanI (|2005l ) and lYoudin &: JohansenI ( 20071 ). The second term on the right 
hand of equation ([50l) also denotes the power due to the fluctuation of dust density. However, the 
dust density fluctuation is produced by the vertical advection under the condition with an initial 
vertical dust density gradient. The fluctuation of dust density is the cause of instability in either 
term of the equation ()50p . The first term arises due to the Lagrangian density fluctuation. On the 
other hand, the cause of the density fluctuation is advection in the second term. The growth rate of 
instability ujj and the phase velocity Vp = ton/k are estimated from the numerical result. Figure [TT] 
shows each term on the left hand of equation ([50j) for TfQx = 10~^, pd{^)/ Pg = 1 and h^jz^ = 0.5. 
Obviously, the energy is gained from the second term associated with Eulerian density fluctuation. 



Because Youdin fc Goodman (j2005l ) and Youdin Sz Johansen ( 2007 ) have assumed that dust 
density is constant, the instability caused by the second term in the equation (f50]l has not been 
seen. In a density stratified layer, if the dust size is small, the instability which obtains the energy 
through the second term in addition to the shear terms occurs. 



m 



Let us see the relatio n between the two-fluid instability and the baroclinic instability stated 
Ishitsu Sekiyal (j2003l ). One fluid approximation assumes infinitesimal friction time rj 0. In 
order for the friction term in equations (llljl and ()13p to be finite, we have Ug — > Ud- Eliminating 
the friction terms by adding equation (fTTI) x{pg/p) and equation (fT3l) x{pd/p) where p = pg + pd, 
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Ud- Linearizing the pressure term in the equation (j5ip yields 
1 d{P + P') 
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p dr p^ dr ^ 
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(52) 



The second term in the right hand is the buoyancy term, which is the cause of the baroclinic 
instability for the one- fluid approximation (jlshitsu fc Sekiyal |2002| ). Note that we assume the gas 
to be incompressible. In the one-fluid model, the fluid of dust and gas mixture is also incompressible 
because V • w = V • li^ = 0. Thus, the dust density perturbation arises not from the compression 
but from the advection; the latter can have non-zero value if 7^ 0. Adding equation (|24p by 
Pg/p and equation ([28|) by Pd/p, we have 



d_ 

dt 



PgUg + 



PdU'd \ 



^ Pg^g ^^'g j_ PdUd du'^ ^ Pg^'g dug _^ Pdw'g du^ 
I p dx p dx p dz I 



dz 



_ PgV'g+PdV'd 
+ ^i'X Z 



p dz 



P 



^SUg 
TfP 



Ud) 



(53) 



It is easily seen that the left hand side and the Coriolis force term in the right hand side become 
those of the linearized equation of equations (|5ip in the limit Ug ^ Ud- We see the frictional term 
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in the right hand of the equation (j53p is written 



' ^'^i^'^^ ^Pi. (54) 



P (Pg + PdY + [pgTfQ.Ky dr 
where we use equations dH), (fT7|) . and (fT9]) . Limiting to the one-fluid leads 



1 _ _ I dP 

K - '^d)p'd -^-^P'd as Tf 0. (55) 

TfPg p'^ or 

This term corresponds to the second term in the right hand of the equation (j52p . However, two-fluid 



instabihty is different from the barochnic instabihty with regards that the barochnic instabihty has 
no axis-symmetric mode. 



5. Conclusions 

We performed the two- fluid of gas and dust, two-dimensional simulations in the dust layer 
of a protoplanetary disk. For TfClx = 10~^, the numerical simulations show the rapid growth 
instability induced mainly by the vertical shear of azimuthal velocity, and additionally the relative 
motion between dust and gas coupled with the dust density fluctuation due to advection if the dust 
density distribution has si gnificant gradient, |^| > /9rf(0)/r/r. The streaming instability stated by 



Youdin &: JohansenI ()2007l ). which is caused by the relative motion of dust and gas coupled with 
the Lagrangian dust density fluctuation, has the small growth rate of the instability if the dust 
size is smaller than several centimeters. On the other hand, the instability powered by the vertical 
shear of the azimuthal velocity, and additionally by the relative velocity of dust and gas coupled 
with dust density fluctuation due to advection shown in this work has the growth rate 0,k even if 
the dust size is small. 

The density fluctuations grows due to the streaming instability if the initial dust density is 
constant, which is accompanied by the concentrations of dust density. However, if the initial dust 
density is not constant, the instability related by the vertical dust density gradient occurs. The 
latter instability diffuses the dust rather than concentrates that. This suggests that the maximum 
density does not always increase from the initial value. 

Additionally, for tj^k = the simulations shows the vertical shear of the radial flow plays 
the important role. After the flow becomes a turbulent state due to the shear instability, the dust 
concentrations are induced by the streaming instability because large relative motion between dust 
and gas is permissible due to loose coupling of dust and gas. 



Chiana (j2008l ) and iBarrancd (j2009l ) have performed the one-fluid, 3D simulations of the shear 
instability and presented that the conditions of the transition to turbulence depends not only on the 
Richardson number but also on the initial perturbations. The dependency of the shear instability 
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on initial conditions is explained as follows. The growth rate of the instability depends on the 
radial and the azimuthal wave numbers k^- and k^ if there is no the radial shear such as the Kepler 
sh ear. The unstabl e regio n in the Fourier space is restricted to a small value of k^- (see figure 1 



m 



Ishitsu &: Sekival (|2003l )). In addition, the shear instability does not have the axis-symmetric 



unstable mode, that is, ks = Q mode. If there is radial shear due to the Kepler motion, the radial 



wave number increases due to shear-stretching as time passes (see equation (43) in llshitsu Sekiya 



(12Q03I)), 

kr{t) = kr{fd) + '^k^tVLK. (56) 

The perturbation can grow only when its wave number passes the unstable region in the Fourier 
space. The flow can transit into turbulence due to the transient amplification if the initial per- 
turbation is large. On the other hand, the flow cannot transit into turbulence for small values of 
initial perturbations. 

However, the instability of two-fluid shown in this work has the axis-symmetric unstable mode. 
As a result, the stabilization caused by the increase of the azimuthal wave number due to the radial 
shear is not effective. We expect that the instability occurs in the radial direction, and then the 
perturbation with small azimuthal wave number grows. 



Chiana (|2008l ) estimates the critical dust surface density as the condition that the shear in- 
stability does not occur if there is not a global turbulence such as MRI. This estimation is derived 
from one-fluid simulations. However, in two-fluid of gas and dust, the instability due to the vertical 
dust density gradient and the relative motion between gas and dust occurs. The flow can transit 



into turbulence even if the disk has the critical surface density estimated by IChiangI (120081 ). Thus, 
the instability induced by the dust density gradient may preclude the planetesimal formation due 
to the gravitational instability. 

In the field of the meteoritics, 1 Myr tim e-lag betwe en the formations of Calcium- Aluminum 
rich inclusions (CAI) and chondrule is known (jscottlboosl ). If the most of planetesimals are formed 
of chondrules, dust aggregates as precursors of chondrules need to be retained in a disk during 1 
Myr. However, the gas friction makes the mm-sized dust in a laminar disk to fall within O.lMyr. 
If the disk is turbulent, some of the dust may avoid falling due to the turbulent diffusion. Thus, 
even though the global turbulence is weak in the dead zone, the turbulence due to the instability 
described in this paper may play the role of avoiding planetesimal formation and floating dust in 
the disk. 

The instability induced by the relative motion between gas and dust should be studied more 
detailedly because this instability have possible important roles in the dust evolution and the 
planetesimal formation in the protoplanetary disk. 



The calculations in this work were partly performed with computers at Astronomical Data 
Analysis Center, National Astronomical Observatory of Japan. This work was supported by Min- 
istry of Education, Culture, Sports, Science and Technology of Japan (MEXT) Grant-in-Aid for 
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Scientific Research on Priority Areas, "Development of Extrasolar Planetary Science" (MEXT- 
16077202). 
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Fig. 1. — Snapshots of dust density at times t^lx = 70 and 85, in the case where TfQx = 1 and 
Pd{0)/pg = 1. Time, length and density are normahzed by Q-k-,W and pg, respectively. 




Fig. 2. — Distributions of dust density, gas and dust velocities of the initial flow in the case where 
TfQx = 10~^ and pd{0)/pg = 1. In the right panel, the curve of the dust velocity is not seen 
because the gas and dust have the almost same azimuthal velocity. 
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Fig. 3. — Snapshots of dust density at times ti^x = 0, 20.0, 24.5 and 29.9, in the case where 
TfQK = 10-3 and Pd{0)/Pg = 1- 




Fig. 4. — The time evolutions of the density weighted averag of the vertical dust r.m.s. defined by 
the equation ([22]) . in the case where tjQk = 10~^(solid) and 1 (dotted) with pd{^)/ pg = 1. 
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Fig. 5. — Distributions of dust density, gas and dust velocities of the initial flow in the case where 
TfQx = 1 and Pd{0)/Pg = 1- 
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6. — Snapshots of dust density at times ti^x = 0, 18.5, 20.5 and 26.5, in the case where 
K = 1 and PdiO)/Pg = 1- 
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Fig. 7. — Each term in the right side of energy equations (|40 p - (j45p . in the case where k'qr = 10.8, 
TfQ.K = 1 and Pd{^)/Pg = 1- 
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Fig. 8. — Each term in the right side of energy equation (j46p . in the case where kr]r = 10.8, 
Tf^K = 1 and Pd{0)/Pg = 1- 



- 23 - 



0.6r 
0.5- 
0.4- 
0.3- 
0.2- 
0.1 - 
0- 



-1 1 1 r 



w 



gas-r 



gas-([) 



gas-z 



i; 



I 



4 -2 



-(dug/dzmw-gu'g*] 

-2Dk5R[v>'/] 



— -idv^dzmw'gv'g*] 

(1/2)Dk5R[«>',*] 

--- -(v^-v,)9?[p>'/]/(z;-p,) 



4-2-10 1 2 3 4 
---- -•3i[{dP'/dz)w'g*]/Pg 

-pM(w'g-w'd)w'g*V(TjPg) 



0.6r 

0.5 
0.4 
0.3 
0.2 
0.1 
0- 



-1 1 1 — 

dust-r 



"1 1 ', 1 1 r- 

dust-(t 



-1 1 1 — 

dust-z 



-2 -1 1 2 3 4 -2 -1 1 2 3 

-Pd(duJdzm[WdU'd*] PjidvJdzmWav'a*] 

-2nK/J^9?[v>'/] -- (l/2)QKyOrf9?[M>'/] 

-P»'a-u'g)u',*]/(TfPg) -p,3i[(v',-v'g)v',*y(TiPg) 



4-2-101 234 

-- -Pj^[(yv'd-^'s>'d*V{TjPg) 



Fig. 9. — Each term in the right side of energy equations (|40 p - (j45p . in the case where krjr = 48.6, 
TfUK = 10-3 and Pd{0)/Pg = 1. 
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Fig. 10. — Each term in the right side of energy equation (j46p . in the case where krjr = 48.6, 
tjQk = 10"^ and pd{0)/pg = 1- 
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Fig. 11. — Each term in the right side of energy equation ()50p . in the case where krjr = 48.6, 
TfVLK = 10-3 and Pd{0)/pg = 1. 
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Table 1. Parameters and results 



run 


Pd(0)/P5 






Nx 


X ^ 


Lx/zd 


Ly/Zd 


smK 




krjr 




const 


1.0 


0.0 


1 


256 


X 256 


TT 


1.0 


5.e-4 


103 


1.9 


0.12 


rlt3 


1.0 


0.5 


10-3 


256 


X 256 


2 TT 


6.0 


5.e-4 


103 


48.6 


1.00 


rltO 


1.0 


0.5 


1 


256 


X 256 


2 TT 


8.0 


2.e-4 


103 


10.8 


1.09 



''The number of Fourier components 



